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Abstract. A new N-body and iiydrodynamical code, called RAMSES, is presented. It has been designed to study 
structure formation in the universe with high spatial resolution. The code is based on Adaptive Mesh Refinement 
(AMR) technique, with a tree based data structure allowing recursiv e grid refinements on a cell-by-cell basis. The 
N-body solver is very similar to the one developed for the ART code (Kravtsov et al. 1997), with minor differences 



in the exact implementation. The hydrodynamical solver is based on a second-order Godunov method, a modern 
shock-capturing scheme known to compute accurately the thermal history of the fluid component. The accuracy 
of the code is carefully estimated using various test cases, from pure gas dynamical tests to cosmological ones. The 
specific refinement strategy used in cosmological simulations is described, and potential spurious effects associated 
to shock waves propagation in the resulting AMR grid are discussed and found to be negligible. Results obtained 
in a large N-body and hydrodynamical simulation of structure formation in a low density ACDM universe are 
finally reported, with 256^ particles and 4.1 x lO'^ cells in the AMR grid, reaching a formal resolution of 8192"^. 
A convergence analysis of different quantities, such as dark matter density power spectrum, gas pressure power 
spectrum and individual haloes temperature profiles, shows that numerical results are converging down to the 
actual resolution limit of the code, and are well reproduced by recent analytical predictions in the framework of 
the halo model. 
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1. Introduction 

Numerical simulations of structure formation in the uni- 
verse are now widely used to understand the highly non- 
linear nature of gravitational clustering. Dark matter is 
believed to be the dominant component in mass of the 
cosmological density field, with only a small fraction (say 
10 %) in baryons. At intermediate scales, such as galaxy 
clusters, dark matter still dominates the total gravita- 
tional mass, but the introduction of a gaseous compo- 
nent appears to be unavoidable, since X-ray or Sunyaev- 
Zeldovich observations of the hot intracluster medium give 
us strong constraints on the structure of galaxy clusters. 
At smaller scales, gas cooling and fluid dynamics play 
a dominant role in the structure of galaxy-size object. 
Although baryons can be described to first order as an 
hydrostatic ionized plasma trapped in dark matter po- 



tential wells, the complexity of hydrodynamical processes 
such as shock heating, atomic radiation cooling and, ulti- 
mately, star formation requires an accurate treatment of 
the baryonic component. 

For a cosmological simulation to be realistic, high mass 
and spatial resolution are needed. While the former is re- 
lated to the initial number of degrees of freedom (usually 
"particles" or "wavelengths") in the computational vol- 
ume, the latter is usually related to the numerical method 
specifically used to compute the particles trajectory. For 
a ri = 1 universe, using a sufficiently large volume of, say, 
100 Mpc h~^ aside, we need at least 256^ particles to de- 
scribe L^, galaxies with 100 particles. In order to resolve 
the internal radial structure of such haloes with at least 
10 resolution elements, we need a spatial resolution of 10 

kpc h~^ or equivalently a dynami cal range of 10^. 

T he Particle-Mesh method ( Hockncy fc Eastwood 



1981 ; Klypin fc Shandarin 1983 ) is perhaps the simplest 



Send offprint requests to: Romain Teyssier 
Correspondence to: Romain.Teyssier@cea.fr 



and fastest N-body algorithm for solving gravitational 
dynamics, but it is limited by computer resources to a 
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dynamical range of 10''. The P3M method (Hockney & 



Eastwc 



od 1981; Efstathiou et al. 1985) can reach a higher 



spatial resolution, by adding a small scale component to 
the PM force, directly computed from the two-body inter- 
actions ("Particle-Particle") between neighboring parti- 
cles. This method suffers however from a dramatic increase 



3- Smooth Particles Hydrodynamics (SPH) can be 
thought as an intermediate solution. This is a particle- 
based method, which follows the Lagrangian evolution 
of the flow, but in which resolution elements are de- 
fined as appropriate averages over 50 neighboring parti- 
cles (Gingold & Monaghan 1977; Evrard 1988; Hernquist 



in CPU time as clustering develops and as short range 
forces become dominant. An improvement of this method 



Katz 1989). This "smoothing kernel" defines the effec- 



was therefore developed for the AP3M code (Couchman 
199l[)T|a hierarchy of recursively refined rectangular grids 
is placed in clustered regions where a local PM solver is 
activated to speed up the PP calculations. Another N- 
body method w hich can achieve high dynamical range is 
the TREE code ( |Barnes fc Hut 1986| ; [Bouchet fc Hernquist 



1988|3vhich properly sort neighboring particles in a recur 
sive tree structure. Long range interactions are computed 
using multipole expansion and low resolution nodes of the 
tree, while short range interactions are computed using a 
PP approach between particles belonging to the same leaf 
of the tree. 

The idea of using Adaptive Mesh Refinement (AMR) 
for N-body methods appears as a natural generalization of 
both AP3M and TREE codes, since they use a hierarchy of 
nested grids to increase the spatial resolution locally. The 
recently developed ART code (Kravtsov et al. 1997) offers, 
in this respect, the first implementation of a grid-based 
high-resolution N-body code, where the mesh is defined 
on a recursively refined spatial tree. ART takes advantage 
of both the speed of a mesh-based Poisson solvers and 
the high-dynamical range and flexibility obtained with 
a tree structure. Since no PP force is considered in the 
ART method, the resolution is not uniform (as opposed 
to AP3M and TREE codes), but proportional to the local 
cell size of the grid. The grid is continuously refined or 
de-refined in the course of the simulation, to ensure that 
the mean number of particles par cell remains roughly 
constant (around 10). This "quasi-Lagrangian" approach 
ensures that two-body relaxation remains unimportant 
( iKnebe et al. 20001 ). 

The gaseous component, baryons, can be described us- 
ing one of several hydrodynamical methods widely used 
today in cosmology. They can be divided into three 
groups: Lagrangian schemes, Eulerian schemes and inter- 
mediate schemes. 

1- Lagrangian schemes or quasi-Lagrangian schemes 
(Gnedin 1995; Pen 1995) are based on a moving mesh 
that closely follows the geometry of the flow for a constant 



numbet of grid pomts. T he grid adapts itself to collapsing 
fluid elements, but sufiers from severe mesh distortion m 



the non linear stage of gravitational clustering ( |Gnedin~fc 
Bertscl linger 199^ ). The coupling with one of the aformen- 
tioned N-body solvers is also non-trivial (ibid). 

2- Eulerian schemes are usually based on uniform 
Cartesian mesh, which make them suitable for a coupling 
with the traditional, low-resolution PM solver (|Cen 1992 ; 



Ryu et] al. 1993|; fleyssier et al. 1998^ [Chieze et al. 199. 
They suffer however from limited dynamical range 



five Eulerian resolution of the method. The SPH method 
offers also the possibility of a straightforward coupling 
to particle-based N-body solver like the AP3M or TREE 
codes. The main drawback of the SPH method is its rel- 
ative poor discontinuity capturing capabilities (one needs 
at least 50 particles per resolution element in order to 
properly describe sharp features, like shock waves or con- 
tact surfaces) and the fact that it relies on the artiflcial 
viscosity method to capture shock waves. 

One of the most promising hydrodynamical method 
at this time is the AMR sch eme, described originall y 
in [Berger fc Ohger (1984| ) and [Berger fc CoUela (1989| ). 



The original AMR method is an Eulerian hydrodynamics 
scheme, with a hierarchy of nested grids covering high- 
resolution regions of the flow. The building blocks of the 
computational grid are therefore rectangular patches of 
various sizes, whose positions and aspects ratio are opti- 
mized with respect to flow geometry, speed and memory 
constraints. Let's call this spatial structure "patch-based 
AMR". An alternative method was proposed by several 
authors (see Khokhlov 199^ , and reference therein) where 
parent cells are refined into children cells, on a cell- by- 
cell basis. As opposed to the original patch-based AMR, 
let us call this last method "tree-based AMR", since the 
natural data structure associated to this scheme is a re- 
cursive tree structure. The resulting grid follows complex 
flow geometry more closely, at the price of a data manage- 
ment which is more complicated than patch-based AMR. 
These two different adaptive mesh structures can be cou- 
pled to any grid-based fluid dynamics scheme. It is worth 
mentioning that modern high-resolution shock capturing 
methods are all grid-based and have number of interest- 
ing features: they are stable up to large Courant num- 
bers, they are striclty conservative for the Euler equations 
and they are able to capture discontinuities within only 
few cells. Among several schemes, higher order Godunov 
methods appear to be more accurate and to be easy to 
generalize in 3 dimensions. 

The original patch-based AMR, based on the Piecewise 
Parabolic Method (PPM: a third-order Godunov scheme), 
was recently adapted to cosmology ( Bryan fc Norman] 
1997; Abel et al. 2000| ). The hydrodynamical scheme was 



coupled to the AP3M N-body solver (without using the 
PP interaction module) . This choice is natural since both 
codes use a set of rectangular patches to cover high- 
resolution regions of the flow. 

In this paper, an alternative solution is explored: cou- 
pling a tree-based AMR hydrodynamical scheme to the N- 



body solver developed for the ART code ( Kravtsov et al 



1997] ). This solution seems indeed more suitable to the 



hierarchical clustering picture where a very complex ge- 
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ometry builds up, with a large number of small clumps 
merging progressively to form large virialized haloes and 
filaments. The number of grids required to cover efficiently 
all these small haloes is so large, that it renders a patch- 
based approach less efficient. 

In this paper, a newly developed N-body and hydrody- 
namical code, called RAMSES, is presented. It is a tree- 
based AMR using the "Fully Threaded Tree" data struc- 
ture of Khokhlov (199S). The N-body solver is largely in- 



spired by the ART code (Kravtsov et al. 1997), with some 
differences in the final implementation. The hydrodynam- 
ical solver is a second-order Godunov scheme for perfect 
gases (also called Piecewise Linear Method or PLM). In 
Section the N-body and hydrodynamical algorithms de- 
veloped for RAMSES are briefly described, with emphasis 
put on the original solutions discovered in the course of 
this work. In Section ||, results obtained by RAMSES for 
standard test cases are presented, demonstrating the ac- 
curacy of the method. Pure hydrodynamical problems are 
considered first, showing that shocks and contact surfaces 
are well captured by the tree-based AMR scheme. 

Great care is taken to demonstrate that refining the 
mesh in shock fronts can be avoided in cosmological con- 
texts. Indeed, potentially spurious effects associated to 
the AMR grid remains low enough to apply the method 
safely using refinements in high density regions only. In 
Section ^, results of a large cosmological simulation using 
256^ particles, with coupled gas and dark matter dynam- 
ics, are reported and compared to various analytical pre- 
dictions. In Section ^, the results presented in this paper 
are summarized, and future projects are discussed. 

2. Numerical methods 

The modules used in RAMSES can be divided into 4 parts: 
the AMR service routines, the Particle Mesh routines, the 
Poisson solver routines and the hydrodynamical routines. 
The dimensionality, noted dim, can be anything among 1, 
2 or 3. 



In RAMSES, time integration can be perfomed in prin- 
ciple for each level independantly. Only two time stepping 
algorithms have been implemented so far: a single time 
step scheme and an adaptive time step scheme. The sin- 
gle time step algorithm consists in integrating the equa- 
tions from t to t + At, with the same time step At for 
all levels. The adaptive time step algorithm, on the other 
hand, is similar to a "W cycle" in the multigrid termi- 
nology. Each level is evolved in time with its own time 
step, determined by a level dependant CFL stability con- 
dition. Gonsequenty, when level £ = is advanced in time 
using one coarse time step, level ^ = 1 is advanced in 
time using two time steps, level £ — 2 using 4 time steps, 
and so on. An additional constraint on these level de- 
pendant time steps comes from synchronization, namely 
AU = Atl^^ + Atj^^. 

Within one time step and for each level, each oper- 
ation is perfomed in the following way: a sub-sample of 
octs is first gathered from the tree. Gathered cells can 
then be modified very efhciently on vector or parallel ar- 
chitectures. Finally, updated quantities are scattered back 
to the tree. When one needs to access neighboring cells (in 
order to compute gradients for example), it is straightfor- 
ward to obtain the neighboring oct addresses from the 
tree, and then to compute the neighboring cell addresses. 

The two main routines used to dynamically modify the 
AMR structure at each time step are now described. 

2.1.1. Building the refinement map 

The first step consists in marking cells for refinement ac- 
cording to user-defined refinement criteria, within the con- 
straint given by a strict refinement rule: any oct in the tree 
structure has to be surrounded by 3'^'™ — 1 neighboring 
parent cells. Thanks to this rule, a smooth transition in 
spatial resolution between levels is enforced, even in the 
diagonal directions. Practically, this step consists in three 
passes through each level, starting from the finer level £max 
down to the coarse grid £ — 0. 



2.1. Adaptive Mesh Refinement 

The fundamental data structu re in RAMSES is called a 
"Fully Threaded Tree" (FTT) ( [Khokhlov 1998| ). Basic ele- 
ments are not single cells, but rather groups of 2"^™ sibling 
cells called octs. Each oct belongs to a given level of re- 
finement labeled £. A regular Gartesian grid, called the 
coarse grid, defines the base of the tree structure (^ = 0). 
In order to access all octs of a given level, octs are sorted 
in a double linked list. Each oct at level £ points to the 
previous and the next oct in the level linked list, but also 
to the parent cell at level ^ — 1, to the 2 x dim neighboring 
parent cells at level £—1 and to the 2'^™ child octs at level 
£ + 1. If a cell has no children, it is called a leaf cell, and 
the pointer to the child oct is set to null. Otherwise, the 
cell is called a split cell. In order to store this particular 
tree structure in memory, one needs therefore 17 integers 
per oct for dim = 3, or equivalently 2.125 integers per cell. 



1. If a split cell contains a children cell that is marked or 
already refined, then mark it for refinement 

2. Mark the 3^^™ - 1 neighboring cells. 

3. If any cell satisfies the user-defined refinement criteria, 
then mark it for refinement. 

One key ingredient still missing in this procedure is the 
so-called "mesh smoothing" . Usually, refinement are acti- 
vated when gradients (or second derivatives) in the flow 
variables exceed a given threshold. The resulting refine- 
ment map tends to be "noisy" , especially in smooth part of 
the flow where gradients fluctuates around the threshold. 
Khokhlov (1998| ) describes a very sophisticated method 
based on a reaction-diffusion operator applied on the re- 
fin ement map. I prefer t o use here the simpler approach 
of Kravtsov et al. (1997 ) in the ART code, where a cu- 
bic buffer is expanded several times around marked cells. 
The number of times one applies the smoothing operator 
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on the refinement map is obviously a free parameter. This 
parameter is noted rtcxpand- In the ART code, this opera- 
tor is appHed twice at each time step. In RAMSES, it is 
usuaUy apphed only once, since, as we see below, bound- 
ary conditions are defined for each level in a slightly more 
sophisticated way than in ART, using buffer regions (see 
§ 2.2.3| ). Therefore, the extra mesh smoothing used in ART 
can be thought as a way of creating the equivalent of the 
buffer regions in RAMSES. 

Note that the exact method implemented here and in 
the ART code leads to a convex structure for the result- 
ing mesh, that is likely to increase the overall stability 
of the algorithm. Note also that only refinement criteria 
are necessary in RAMSES: no de-refinement criteria need 
to be specified by the user. This is an important differ- 



A collisionless N-body system is described by the 
Vlasov-Poisson equations, which, in terms of particles (la- 
beled by "p")j reads 



ence compared to other approaches (Kravtsov et al. 1997; 
Khokhlov 19"9^ ). 



2.1.2. Modifying the tree structure 

The next step consists in splitting or destroying children 
cells according to the refinement map. RAMSES performs 
two passes through each level, starting from the coarse 
grid £ = 0, up to the finer grid £max: 

1. If a leaf cell is marked for refinement, then create its 
child Oct. 

2. If a split cell is not marked for refinement, then destroy 
its child Oct. 

Creating or destroying a child oct is a time-consuming 
step, since it implies reorganizing the tree structure. 
Thanks to the double linked list associated to the FTT 
tree structure, this is done very efficiently by first discon- 
necting the child oct from the list, and then reconnect- 
ing the list in between the previous and next octs. Note 
however that this operation can not be vectorized. It is 
important to stress that this operation is applied at each 
time step, but for a very small number of octs. In other 
word, at each time step, the mesh structure is not rebuilt 
from scratch, but it is slightly modified, in order to follow 
the evolution of the flow. Since the refinement map has 
been carefully built during the last step, the refinement 
rule should be satisfied by construction. This is however 
not the case if one uses the adaptive time step method 
described in Section 2.3.2. In this final check is 



performed before splitting leaf cells. If the refinement rule 
is about to be violated, leaf cells are not refined. 



2.2. N-body solver 

The N-body scheme used in RAMSES is similar in many 



aspects to the ART code of Kravtsov et al. (1997). Since 
the ART code was not publicily available at the time this 
work was initiated, a new code had to be implemented 
from scratch. I briefly recall here the main ingredients of 
the method, outlining the differences between the two im- 
plementations. 



— - — v„ and — - — —Vxip where Aj-t/i = AirGp (1) 
at at 



Grid-based N-body schemes, such as the standard PM, 
are usually decomposed in the following steps: 

1. Compute the mass density p on the mesh using a 
"Cloud-In-Cell" (CIC) interpolation scheme. 

2. Solve for the potential (p on the mesh using the Poisson 
equation. 

3. Compute the acceleration on the mesh using a stan- 
dard finite-difference approximation of the gradient. 

4. Compute each particle acceleration using an inverse 
CIC interpolation scheme. 

5. Update each particle velocity according to its acceler- 
ation 

6. Update each particle position according to its velocity. 

The specific constraints of a tree-based AMR N-body 
solver are now discussed in more details. 

2.2.1. The particle linked list 

Since we are dealing with an AMR grid, we need to know 
which particle is interacting with a given cell. This is done 
thanks to a particle linked list. Particles belong to a given 
oct, if their position fits exactly into the oct boundaries. 
All particles belonging to the same oct are linked together. 
In order to build this linked list, we have to store the 
position of each oct in the tree structure. Moreover, each 
oct needs to have access to the address of the first particle 
in the list and to the total number of particles contained 
in its boundaries. We need therefore to store these two 
new integers in the FTT tree structure. 

The particle linked list is built in a way similar to the 
TREE code: particles are first divided among the octs sit- 
ting at the coarse level £ = 0. Each individual linked list 
is then recursively divided among the children octs, up 
to the finer level £ = £max- Going from level £ to level 
i -\- 1 implies removing from the linked list particles sit- 
ting within split cell boundaries. Going from level £ -\- 1 
to level £ implies reconnecting the children linked lists to 
the parent one. In the adaptive time step case, in order to 
avoid rebuilding the whole tree from the coarse level, par- 
ticle positions are checked against parent octs boundaries 
and, if necessary, are passed to neighboring octs using the 
information stored in the FTT tree. 

2.2.2. Computing the density field 

The density field is computed using the CIC interpola- 
tion scheme ( Hockney fc Eastwood 1981 ). For each level, 
particles sitting inside level £ boundaries are first consid- 
ered. This can be done using the level £ particle linked 
list. Particles sitting outside the current level, but whose 
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clouds intersect the corresponding volume are then taken 
into account. This is done by examining particles sitting 
inside neighboring parent cells at level £ — 1. Note that 
in this case the size of the overlapping cloud is the one 
of level £ particles. In this way, for a given set of particle 
positions, the resulting density field at level £ is exactly 
the same as that of a regular Cartesian mesh of equivalent 
spatial resolution. 

2.2.3. Solving the Poisson equation 

Several methods are described in the literature to solve 
for the Poisson equation in the adaptive grids framework 



(Couchman 1991; Jessop ct al. 1994; Kravtsov et al. 1997; 



Almgren et al. 1998; Truelove et al. 1998). In RAMSES 



as in the ART code, the Poisson equation is solved using a 



"one-way interface" scheme (Jessop et al. 1994; Kravtsov 



et al. 1997): the coarse grid solution never "sees" the ef- 



fect of the fine grids. The resulting accuracy is the same 
as if the coarse grid were alone. Boundary conditions are 
passed only from the coarse grid to the fine grid by a lin- 
ear interpolation. For each AMR level, the solution should 
therefore be close to the one obtained with a Cartesian 
mesh of equivalent spatial resolution, but it can not be 
better in any way. A better accuracy would be obtained 
using a two- way interface scheme, as the one described for 



example in Truelove et al. (19981) . Such a sophisticated 



improvement of the Poisson solver is left for future work. 

The Poisson equation at the coarse level is solved us- 
i ng standard Fast Fourier Transform (FFT) technique 
( Hockney fc Eastwood 1981 ), with a Green function cor- 
responding to Fourier transform of the 2dim -I- 1-points 
finite difference approximation of the Laplacian. For fine 
levels {£ > 0), the potential is found vising a relaxation 
method similar to the one developed for the ART code: the 
Poisson equation is solved using the 2dim + 1-points finite 
difference approximation of the Laplacian, with Dirichlet 
boundary conditions. In RAMSES, boundary conditions 
are defined in a temporary buffer region surrounding the 
level domain, where the potential is computed from level 
£ — 1 through a linear reconstruction. 

Using these specific boundary conditions, the potential 
can be computed using any efficient relaxation method. In 
RAMSES, the Gauss-Seidel (GS) method with Red-Black 



Ordering and Successive Over Relaxation (Press et 



al. 



1992) is used. In two dimensions, for unit mesh spacing. 



the basic GS writes as 



1 

4 



,-i)-4a. (2) 



This iteration is applied first to update the potential for 
"black" cells defined by i odd and j odd or i even and j 
even, and then to update the potential for "red" cells de- 
fined by i odd and j even or i even and j odd. Finally, 
the result is "over-corrected" using the so-called over- 
relaxation parameter uj 



The speed of the algorithm relies on the correct choice for 
both UJ and the initial guess (fi^j. For a regular N x N 
Cartesian mesh, the optimal over-relaxation parameter is 
known to be (Press et al. 1992) 



UJ 



1 



(4) 



with 1 <uj <2 



(3) 



where a = 1 for Dirichlet boundary conditions and a = 2 
for periodic boundary conditions. For an irregular AMR 
grid, the situation is more complicated, since the compu- 
tational volume is covered by irregular mesh patches. The 
over-relaxation parameter has to be found empirically. An 
interesting way of determining the optimal value for uj is 
to estimate the average size < L > of these patches, and 
to use it in formula (^) in place of N. It was found to work 
very well in practice. 

The initial guess is obtained from the coarser level £ — 1 
through a linear reconstruction of the potential. In this 
way, the solution at large scale is correctly captured at the 
very beginning of the relaxation process. Only the shortest 
wavelengths need to be further corrected. 

A question that arises naturally is: when do we reach 
convergence ? Since our Poisson solver is coupled to a N- 
body system, errors due to the CIC interpolation scheme 
are dominant in the force calculation. As soon as the resid- 
uals are smaller than the CIC induced errors, further iter- 
ations are unnecessary. For cosmological simulations, this 
is obtained by specifying that the 2-norm of the residual 
has to be reduced by a factor of at least 10'^. 

Let us consider a 128'^ coarse grid, completely refined 
in a 256'^ underlying fine grid. Solving first the Poisson 
equation on the coarse level using FFT, the solution is 
injected to the fine grid as a first guess. In this particular 
example, the optimal over-relaxation parameter is w ~ 1.9 
(using Eq. ^ with a = 2) and 60 iterations are needed to 
damped the errors sufficiently. 

Let us now consider a more practical example, in 
which a typical AMR grid is obtained from a cosmologi- 
cal simulation. In this case, the average AMR patch size 
was empirically found to be roughly < L >~ 20 cells. 
Equation (^ with a = 1 gives uj ~ 1.7. 20 iterations only 
are needed for the iterative solver to converge sufficiently. 
Note that for the ART code, the optimal value was found 
to be w = 1.25 (Kravtsov et al. 1997), using however a 
different approach to set up boundary conditions. 



2.2.4. Computing the acceleration 

Using the potential, the acceleration is computed on the 
mesh using the 5-points finite difference approximation of 
the gradient. As for the potential, the acceleration is cell- 
centered and the g radient stencil is symmetric al in order 
to avoid self- forces ( Hockney fc Eastwood 1981 ). Buffer re- 
gions defined during the previous step are used here again 
to give correct boundary conditions. The acceleration is in- 
terpolated back to the particles of the current level, using 
an inverse CIC scheme. Only particles from the linked list 
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whose cloud is entirely included within the level bound- 
ary are concerned. For particles belonging to level £, but 
whose cloud lies partially outside the level volume, the ac- 
celeration is interpolated from the mesh of level £ — 1. This 
is the same for the ART code: "In this way, particles are 
driven by the coarse force until they move sufficiently far 
into the finer mesh" (Kravtsov et al. 1997). 



2.2.5. Time integration 

One requirement in a coupled N-body and hydrodynami- 
cal code is the possibility to deal with variable time steps. 
The stability conditions for the time step is indeed given 
by the Courant Friedrich Levy (CFL) conditio n, which 
can vary in time. T he standard leapfrog scheme ( Hockney 
fc Eastwood 1981 ), though accurate, does not offer this 
possibility. In RAMSES, a second-order midpoint scheme 
has been implemented, which reduces exactly to the sec- 
ond order leapfrog scheme for constant time steps. Since 
the acceleration — V(/)" is known at time t" from particle 
positions Xp , positions and velocities are updated first by 
a predictor step 



and then by a corrector step 



n+l ^ n+1/2 



(5) 
(6) 



(7) 



In this last equation, the acceleration at time is 
needed. In order to avoid an extra call to the Poisson 
solver, this last operation is postponed to the next time 
step. The new velocity is computed as soon as the new 
potential is obtained. In RAMSES, it is possible to have 
either a single time step for all particles, or individual time 
steps for each level. In the latter case, when a particle exits 
level £ with time step Ati, the corrector step is applied at 
level £—1, using Atg in place of Ai^-i. Therefore, the "past 
history" of all particles has to be known in order to apply 
correctly the corrector step. This is done in RAMSES by 
introducing one extra integer per particle indicating its 
current level. This particle "color" is eventually modified 
at the end of the corrector step. 

Usually, the time step evolution is smooth, making our 
integration scheme second-order in time. However, if one 
uses the adaptive time step scheme instead of the more ac- 
curate (but time consuming) single time step scheme, the 
time step changes abruptly by a factor of two for particles 
crossing a refinement boundary. Only first order accuracy 
is retained along those particle trajectories. This loss of 
accurac y has been analyzed in realistic cosmological con- 
ditions ([Kravtsov fc Klypin 1999i [Yahagi fc Yoshii 200lD 



and turns out to have a small effect on the particle distri- 
bution, when compared to the single time step case. 



2.3. Hydrodynamical Solver 

In RAMSES, the Euler equations are solved in their con- 
servative form: 



dt 



| + V.(pu)^0 
(pu) + V • (pu (g) u) + Vp = -pS/(j) 



dt 



(pe) -I- V • [pu (e -I- p/ p)] = -pu • V(?!) 



(8) 
(9) 

(10) 



where p is the mass density, u is the fluid velocity, e is the 
specific total energy, and p is the thermal pressure, with 



P = (7 - 



(11) 



Note that the energy equation (Eq. |T^) is conservative 
for the total fluid energy, if one ignores the source terms 
due to gravity. This property is one of the main advan- 
tages of solving the Euler equations in conservative form: 
no energy sink due to numerical errors can alter the flow 
dynamics. Gravity is included in the system of equation 
as a non stiff source term. In this case, the system is not 
explicitly conservative and the total energy (potential -|- 



kinetic) is conserved at the percent level (see section 4.3) 



Let [/" denote a numerical approximation to the cell- 
averaged value of (p, pu, pe) at time and for cell i. The 
numerical discretization of the Euler equations with grav- 
itational source terms writes: 



n+l 



UT 



pn+1/2 _ pn+1/2 



i+1/2 



1/2 _ ^n+1/2 



At 



Ax 



(12) 



The time centered fluxes F. 



1+1/2 
+1/2 



across cell interfaces are 



computed using a second-order Godunov method (also 
known as Pieceweise Linear Method), with or without di- 
rectional splitting (according to the user's choice), while 
gravitational source terms are included using a time cen- 
tered, fractional step approach: 



ipu)-y<f>f + {pu)-+'vcj^^i 



A general description of Godunov a nd fractional step 
methods can be found in fPoro (19971). The present im- 



plemcntation is based on the work of [CoUcla (19901) and 
zman (1994[ ). For sake of brevity, only its basic fea- 



tures are recalled here. 

2.3.1. Single grid Godunov solver 

In this section, I describe the basic hydrodynamical 
scheme used in RAMSES to solve equations (p[]lO|) at a 
given level. It is assumed that proper boundary conditions 
have been provided: the hydrodynamical scheme requires 
2 ghost zones in each side and in each direction, even in 
the diagonal directions. Since in RAMSES the Euler equa- 
tions are solved on octs of 2"^™ cells each, 3*^"° — 1 similar 
neighboring octs are required to define proper boundary 
conditions. The basic stencil of the PLM scheme therefore 
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contains 6'^™ cells. This is not the case for PPM ( goUela 



& Woodward 1984) for which 4 ghost zones are required in 
each side and in each direction. Since the AMR structure 
in RAMSES is based on octs (2^^™ cells), PPM would be 
to expensive to implement in many aspects. One solution 
would be to modify the basic tree element and increase the 
number of cells per oct from 2'^™ cells to 4^^™ cells. The 
resulting AMR structure would however loose part of its 
flexibility to adapt itself to complex flow geometry. The 
FLASH code (Fryxell et al. 2000) is an example of such 
an implementation, using the PPM scheme in a similar 
recursive tree structure, with however S'^™ cells per basic 
tree element. 

For a given time step, we need to compute second- 
order, time-centered fluxes at cell interfaces. This is done 
in RAMSES using a Riemann solver, with left and right 
states obtained by a characteristics tracing step. A stan- 
dard characteristic analysis is done first, by Taylor ex- 
panding the wave equations to second order and projecting 
out the waves that cannot reach the interface within the 
time step. These states are then adjusted to account for 
the gravitational field. If the chosen scheme is not direc- 
tionnaly splitted, transverse derivative terms are finall y 
added to account for transverse fluxes ( Saltzman 1994 ). 
The slopes that enter into the Taylor expansion are com- 
puted using the Min-Mod limiter to ensure the monotonic- 
ity of the solution. 

The Riemann solver used to compute the Godunov 
states is "almost exact" , in the sense that a correct pres- 
sure at the contact discontinuity is obtained iteratively 
(typically, for strong shocks 10 Newton-Raphson itera- 
tions are required for single-precision accuracy of 10^^). 
The only approximation relies in the assumption that the 
rarefaction wave has a linear proflle. In the final step, 
fluxes of the conserved variables are computed using these 
Godunov states. The outputs of the single grid algorithm 
are therefore fluxes across cell interfaces. 

Practically, this single grid module is applied to a large 
vector of stencils of 6'^™ cells each. For a large Cartesian 
grid of A^'^™ cells, the CPU time overhead associated to 
this solution is rather large. Since the main time consum- 
ing part is the Riemann solver, the estimated CPU time 
overhead was found to be roughly 50 %, 100 % and 200 % 
for dim=l, 2 and 3 respectively. Since in any useful AMR 
calculation, the mesh structure is not a regular Cartesian 
grid, the actual overhead is much lower, although difficult 
to estimate in practice. Moreover, this solution is much 
easier to implement than any potentially faster alterna- 
tive one can think of, and easy to optimize on vector and 
parallel supercomputers. 

2.3.2. AMR implementation 

This section describes how the solution is advanced in time 
within the present AMR methodology. Note that this pro- 
cedure is recursive with respect to level i (step 3). 



1. Generate new refinements at level l+lhy conservative 
interpolation of level I variables. 

2. Compute the new time step At^ using the GEL 
Courant condition and the constraint At^ < A<£_i. 

3. Advance the solution in time for level ^ -I- 1, once in 
the single time step case, or twice for the adaptive time 
step case. 

4. Modify the time step At^ according to the synchroni- 
sation constraint A<£ = At^+i for the single time step 
case or Ai^ ~ A<)^j^ + Ai^_|_j^ for the adaptive time step 
case. 

5. Compute boundary conditions in a temporary buffer 
by conservative interpolation of level i — 1 variables. 

6. Compute fiuxes using the single grid Godunov solver. 

7. Replace the fluxes at coarse-fine interface by averaging 
the fluxes computed at level £ + 1. 

8. For leaf cells, update variables using these fluxes. 

9. For split cells, update variables by averaging down the 
updated variables of level i + 1. 

10. Build the new reflnement map. 



In RAMSES, boundary conditions are supplied to flne lev- 
els by a conservative linear reconstruction of coarse cell 
values (step 5). The actual interpolatio n scheme is a 3D 
generalization of the Min-Mod limiter ( De Zeeuw 1993 ). 
The coarse solution is assumed to remain constant in time 
during the advance of the fine solution. For flne cells at 
coarse-flne boundaries and for the adaptive time step case 
only, the accuracy reduces thus from second to first or- 
der in time, but the global solution remains second order 
(iKhokhlov 1998|). 



2.4. Time Step Control 

The time step is determined for each level independently, 
using standard stability constraints for both N-body and 
hydrodynamical solvers. 

The first constraint comes from the gravitational evo- 
lution of the coupled N-body and hydrodynamical sys- 
tem, imposing that At^ should be smaller than a fraction 
Ci < 1 of the minimum free-fall time 



At{ = Ci X min(tff) 



(14) 



An additional constraint comes from particle dynamics 
within the AMR grid, imposing that particles move by 
only a fraction C2 < 1 of the local cell size. 



Aij = C2 X Ax / max{vp) 



(15) 



A third constraint is imposed on the time step by speci- 
fying that the expansion factor Cexp should not vary more 
than C3 ~ 10 % over one time step. This constraint is ac- 
tive only at early times, during the linear regime of grav- 
itational clustering. 



A4 = C3 X flexp/ 



(16) 
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The last constraint is imposed by the Courant Friedrich 
Levy stabihty condition, which states that the time step 
should be smaller than 

A4 = cflxA.7max(K|+c,K|+c,K|+c) (17) 

where cfl < 1 is the Courant factor. In the coupled N-body 
and hydrodynamics case, the actual time step is equal to 
min(Aif , A4, A4, A4). 

2.5. Refinement Strategy 

The refinement strategy is the key issue for any AMR 
calculation. Bearing in mind that the overhead associated 
to the AM R sch eme can be as large as a factor of 2 to 3 
(see section 2.3) compared to the corresponding uniform 
grid algorithm, the maximum fraction of the grid that 
can be refined lies in between 30 % to 50 %. One should 
therefore design a refinement strategy that allows for an 
accurate treatment of the underlying physical problem, 
but minimizes also the fraction of the volume to be refined. 

For the N-body solver, the refinement strategy is 
based on the so-called "quasi-Lagrangian" approach. As 
in Kravtsov et al. (1997), the idea is to obtain a constant 



number of particles per cell. In this way, two-body relax- 
ation effects can be minimized, as well as the Poisson noise 
due to particle discreteness effects. The latter effect can 
be damaging when coupling the N-body code to the hy- 
drodynamics solver. The "quasi-Lagrangian" approach is 
implemented by refining cells at level £ if the dark mat- 
ter density exceeds a level dependent density threshold, 
defined as 



pe = Mc X {Ax") 



e\-diia 



(18) 



where Mc is the maximum mass (or number of particles) 
per cell. For pure N-body simulations, Mc is usually cho- 
sen around 5-10 particles ( Kravtsov ct al. 1997| ), which 
gives a few particles per cell on average. For gas dynamics 
simulations, Mc should be chosen around 40-80 particles, 
in order to lower enough the Poisson noise. In this case, 
we obtain indeed more than 10 particles per cell on av- 
erage. Note also that since for gas dynamics simulations, 
the total memory is dominated by the storage associated 
to the fiuid variables, the number of particles par cell can 
be chosen much higher than for pure dark matter simula- 
tions. 

As for the N-body solver, a "quasi-Lagrangian" ap- 
proach can be implemented for the gas component, using 
level dependent density thresholds defined by 



Pi = Mb X (Ax^ 



>-dii 



(19) 



In order to follow the same Lagrangian evolution than the 
dark matter component, the typical baryonic mass per cell 
Alb can be derived as 



Mh = Mc 



For pure gas dynamics applications, other refinement cri- 
teria can be used (see Khokhlov 199? , for more exam- 
ples). In RAMSES, only refinement criteria based on gra- 
dients of the fiow variables have been implemented: for 
each cell i and for any relevant flow variable q (pressure, 
density, Mach number...), its gradient is computed using 
the 2 X dim neighboring cells. If this gradient, times the 
local mesh spacing, exceeds a fraction of the central cell 
variable 



Aa 



(21) 



then cell i is refined. The parameter Cq is a free parameter 
that need to be specified by the user. A similar criterion 
based on second derivatives of the flow variables has also 
been implemented. 

The last refinement criterion implemented in RAMSES 
is purely spatial: for each level, refinements are not allowed 
outside a sphere centered on the box center. This last cri- 
terion allows the user to refine the computational mesh 
only in the center of the box, in order to follow properly 
the formation of a single structure, without spending to 
much resources in refining also the surrounding large scale 
field. The radius of this spherical region, noted Rg, can be 
specified for each level independently. 

2.6. Cosmological Settings 

RAMSES can be used for standard fiuid dynamics or N- 
body problems, with periodic, reflecting, inflow or outflow 
boundary conditions. For the present paper, RAMSES 
is however used in the cosmological context. The N- 
body solver and the hydrodynamics solver are both im- 
plemented using "conformal time" as the time variable. 
This allows a straightforward implementation of comov- 
ing coordinates, with minor changes to the original equa- 
tions. The details of these so-called "super-comoving co- 
ordinates" can be found in Martel & Shapiro (1998) and 



references therein. The idea is to perform the following 
change of variables 



di = Ho- 



dt 



and X ~ 



1 X 
a L 



^mPc 



HnL 



and P 



P 



(22) 



(23) 



(24) 



where is the Hubble constant, f7,„ is the matter density 
parameter, L is the box size and pc is the critical density. 
In the specific case 7 = 5/3, equations (|l|) and (p[]lO|) re- 
main unchanged, at the exception of the Poisson equation 
which now reads 



A, 



-aVlm{p - 1) 



(25) 



If 5/3 a single additional source term must be in- 

(20) 

eluded in the right-hand side of the energy conservation 
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equation (Eq. ^0|). These "super-comoving coordinates" 
simplify greatly the introduction of comoving variables in 
the equations. 

3. Tests of the code 

In this section, I present tests of increasing complexity for 
both the N-body solver and the hydrodynamical solver. 
These tests are also useful to choose the correct parame- 
ters for realistic cosmological applications described in the 
last section. 

3.1. Acceleration around a Point Mass 

Particles of unit mass are placed randomly in the compu- 
tational box, whose coarse grid is defined by nx = ny = 
riz = 32. Test particles are then dropped randomly in or- 
der to sample the acceleration around each massive par- 
ticle. The AMR grid is built around each central particle. 
For that purpose, refinement density thresholds were set 
io pp = Q for each level. An increasing number of refine- 
ment levels was used, from ^max = to ^max — 6, the 
latter case corresponding to a formal resolution of 2048'^. 
Mesh smoothing was performed with rtcxpand = 1- 

Figure ^ shows the resulting radial and tangential ac- 
celerations, divided by the true force. The tangen- 
tial acceleration gives here an indication of the level of 
force anisotropy and accuracy. Note that the acceleration 
due to the ghost images of the massive particle (periodic 
boundary conditions) was substracted from the computed 
acceleration (using the Ewald summation method). For 
comparison, the acceleration of an homogeneous sphere 
(with radius equal to the cell size of the maximum refine- 
ment level) is also plotted in Figure |l| as a solid line : 
the AMR acceleration appears to provide a slightly lower 
spatial resolution (rouhly 1.5 cell size). At lower radius, 
the force smoothly goes to zero, exactly as for a PM code 
of equivalent dynamical range. At lower radius, the force 
anisotropy is also the same as for a PM code of equivalent 
dynamical range. On the other hand, at higher radius, the 
force anisotropy remains close to 1 %. Contrary to a single 
grid PM solver, the force error does not decrease mono- 
tonically as radius increases. Here, the error level remains 
roughly constant (at the percent level), since the spatial 
resolution also decreases as radius increases. In fact, the 
AMR force on a given particle corresponds to a single grid 
PM force whose cell size is equal to that of the particle's 
current level. As one goes from one level to the next, dis- 
continuities in the force remain also at the percent level. 

3.2. Acceleration of Particles in a ACDM Simulation 

In order to assess the quality of the gravitational accel- 
eration computed by RAMSES in a cosmological situa- 
tion, we consider now a set of 64'^ particles obtained in 
a ACDM simulation, at redshift z = 0. In this way, we 
are able to quantify the force errors in a typical hier- 
archical clustering configuration, with the corresponding 



mesh refinements structure. The coarse grid was defined 
by Tlx — riy — Uz — 32 and each particle was assigned a 
mass TUp — 1/8. The adaptive mesh was built using re- 
finement density thresholds — 5 x 8^ for each level £. 
Each cell is therefore refined if it contains more than 40 
particles, with a roughly constant number of particles per 
cell after each refinement (between 5 and 40) . Mesh struc- 
tures associated to this particle distribution are shown in 
the last section of the paper (Fig. pO| ). 

For each level of refinement, the AMR force is then 
compared to the PM force of equivalent spatial resolu- 
tion (see Fig. ||). For particles sitting at the coarse level 
{£ = 0), the force is by construction exactly equal to the 
PM force with 32^^ cells (results not shown in the figure). 
For levels ^ = 1, ^ = 2, f = 3 and £ = 4, the AMR force 
is compared to the PM force with respectively 64'^, 128^, 
256'^ and 512^ cells. In Figure ||, each panel shows the dif- 
ference between the AMR force and the PM force for each 
level. The number of particles sitting at each level is indi- 
cated in the upper-left part of each panel. The mean force 
error and the standard deviation is indicated in the lower- 
left part of each panel. Although the error distribution is 
strongly non Gaussian, its typical magnitude remains at 
the percent level in all cases. Note that errors are larger 
for forces of intermediate and small values, indicating that 
those particles might be sensitive to the boundary condi- 
tions (tidal field) imposed on the level boundaries (this is 
the main source of inaccuracy in the N-body scheme). On 
the other hand, for particles with strong acceleration, the 
AMR force is almost indistinguishable from the PM force 
of equivalent resolution. 

3.3. Shock Tube 

The initial conditions are defined by a left state given by 
PL = I, ul = and Pl — I and by a right state given 
by PR = 0.125, UR = and Pr = 0.1 for a 7 = 1.4 
fiuid. This test (also called Sod's test) is interesting be- 
cause it captures all essential features of one dimensional 
hydrodynamical flows, namely a shock wave, a contact dis- 
continuity and a rarefaction wave. While the latter wave 
remains continuous, the 2 former features are discontinu- 
ous. Modern shocks capturing methods like the one used 
in RAMSES spread shock fronts over 2 to 3 zones. Contact 
discontinuities are usually more difficult to capture (say 6 
to 10 cells), and the spreading usually increases with the 
number of time step. This numerical smoothing is respon- 
sible for the dissipation of the scheme. AMR technique al- 
lows one to increase the spatial resolution around the dis- 
continuities and therefore to minimize the numerical dissi- 
pation. In the present application, the refinement criteria 
are based on p ressure, density and Mach number gradients 
(see Sect, p^ ), with parameter Cp = Cp — Cm = 0.01. 
The maximum number of refinements was set to i'max = 6, 
fo r a co arse level mesh size rix = 64. Mesh smoothing (see 
§ |2.1.l|) is performed using noxpand = 1- Note that the 
refined mesh is built before the beginning of the Simula- 
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Fig. 1. Acceleration of massless test particles dropped randomly around single massive particles, whose positions 
are also chosen randomly in the box. The coarse grid has 32'^ cells. The number of refinement levels is progressively 
increased from to 6, with increased spatial resolution around the massive particles. The radial AMR acceleration 
divided by the true acceleration is shown as light grey dots (versus radius in units of coarse cell size). The same ratio 
for the tangential AMR acceleration is shown as dark grey dots. The force corresponding to an homogeneous sphere 
with radius equal to the smallest cell length is also plotted for comparison (solid line). 

Fig. 2. Particle positions obtained in a ACDM simulation are considered in this test. Each panel shows the force error 
for particles sitting at different levels of refinement. The error is defined as the difference between the AMR force and 
the force computed by a PM code of equivalent spatial resolution. In each panel, the average and the variance of the 
error are also shown. 

Fig. 3. Shock Tube Test: Density, velocity, pressure and refinement level as a function of position at time t — 0.245. 
Numerical results are shown as squares, and compared to the analytical solutions (solid lines). See text for details. 



tion. The time step is controlled by a Courant number 
cfl = 0.8. Results are shown at time t = 0.245 and com- 
pared to the analytical solution in Figure |[ The shock 
front and the contact surface are refined up to the max- 
imum refinement level: the formal resolution is therefore 
4096. The total number of cells (counting both split and 
leaf cells) is only 560, or 14 % of the corresponding uni- 
form mesh size. 69 time steps were necessary at the coarse 
level, while 4416 time steps were necessary at level £ = 6. 
It is worth mentioning that pressure and velocity remain 
remarkably uniform across the contact discontinuity, and 
no side effects due to the presence of discrete refinement 
ratio are noticeable. 



3.4. Planar Sedov Blast Wave 

The last test, though interesting and complete, is not a 
very stringent one, since it involves a rather weak shock. 
In order to test the ability of RAMSES to handle strong 
shocks (a common feature in cosmology), let us consider 
the planar Sedov problem: the computational domain is 
filled with a 7 — 1.4 fluid with po — I, uq — and 
Pq = 10~^. A total (internal) energy Eq = 1/2 is deposited 
in the first cell only at x = 0+. Note that here again the 
refined mesh is built before the beginning of the simula- 
tion. Reflexive boundary conditions are considered. The 
grid is defined by — 32 with 6 levels of refinement. The 
only refinement criterion used here is based on pressure 
gradients, with Cp = 0.1. Mesh smoothing is guaranteed 
by 'T-oxpand = 1- The Courant number is set to cfl = 0.8. 
Very rapidly, a self-similar flow builds up, following the 
analytical solution described in Sedov (1993). Simulation 
results are shown for different output times and compared 
to the analytical solutions. Note that the shock front prop- 
agates exactly at the correct speed. The numerical solu- 
tion closely matches the analytical one, without any visible 
post-shock oscillations. 239 time steps only were necessary 
at the coarse level, but 15296 time steps were necessary 
at the finest refinement level. The total number of cells 
(including split cells) in the adaptive mesh structure is 
only 130, to be compared with 2048 cells for the uniform 



Fig. 5. Strong Shock Passing through a Coarse-Fine 
Interface with Mach number M = 5000 and 7 = 5/3, com- 
puted with cfl — 0.5. The upper plot shows the reference 
case with a 256 cells coarse grid, uniformly refined up to 
level i — 1 (without any coarse-fine boundary). The mid- 
dle figure shows the case where the shock goes through a 
fine-to-coarse boundary (located around x ~ 234), while 
the bottom figure shows the case where the shock goes to a 
coarse-to-fine boundary (located around the same place). 
In the latter case, perturbations of the order of 5 % are 
generated in the post-shock flow. 



grid of equivalent spatial resolution (4.3 %). Due to the 
refinement criterion used here, the adaptive mesh mainly 
concentrates the computational effort around the shock 
front. In one dimension, as it is the case here, disconti- 
nuities like shocks are quite inexpensive to deal with: if 
one adds one level of refinement, the total number of cells 
increases by a constant (and small) amount. For two- and 
three-dimensional calculations, the situation is much more 
demanding: since shocks and contacts discontinuities are 
surface waves, increasing the resolution by a factor of 2 
corresponds to increasing the total number of cells by a 
factor of 2 for dim = 2 and 4 for dim = 3. Therefore, we 
have to face the possibility of stopping at some level the 
refinement hierarchy and investigate what happens to the 
numerical solution in this case. 

3.5. Strong Shock Passing through a Coarse-Fine 
Interface 

It is well known that in any AMR calculations, disconti- 
nuities in the flow (like the one discussed in the previous 
sections) must be refined up to the maximum level, in 



order to obtain accurate results (Berger & CoUcla 1989; 
[Khokhlov 199S ) . Unfortunately, it is not always possible to 



satisfy this rule because of memory limitations, even on 
modern computers. One has therefore to consider cases 
for which discontinuities leave or enter regions of different 
spatial resolution. The situation is especially sensitive for 
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Fig. 4. Planar Sedov Blast Wave Test: Density, velocity, pressure and refinement level as a function of position for 
times t =0.05, 0.1, 0.15, 0.2, 0.25 and 0.3. Numerical results are shown as squares, and compared to the analytical 
solutions (solid lines). See text for details. 



contact surfaces, since as soon as the code spreads them 
over, say, 6 cells, no matter how much one refines them af- 
terwards, they will preserve their original thickness. Shock 
waves, however, have a self-steepening mechanism that al- 
lows them to adapt to the local resolution and restore their 
sharp profile over 2 to 3 cells only. The price to pay for this 
interesting property is the appearance of post-shock oscil- 
lations after the front has entered a high-resolution region. 
To illustrate this, Khokhlov (199^ ) proposed a simple test 
based on the propagation of a strong shock wave across a 
coarse-fine interface. Khokhlov's test is reproduced here 
using the following parameters: the shock Mach number 
is set to M = 5000 with 7 — 5/3. The base grid resolu- 
tion is set to Tlx = 256 and the Courant number is set to 
cfl = 0.5. 

The following three cases have been considered: 

1. the whole computational domain is refined up to £ = 1. 

2. the computational domain is refined up to ^ = 1 only 
left to x = 234. 

3. the computational domain is refined up to ^ = 2 right 
to 2: = 234 and up to = 1 otherwise. 

The resulting density profiles are shown in Figure |[ While 
in the 2 former cases, the density profiles show no vis- 
ible oscillations, the latter case does show small oscilla- 
tions of the order of 5 %. This is a direct consequence of 
the abrupt change of spatial resolution between the 2 lev- 



the total number of cells necessary to cover the shock front 
scales with spatial resolution as 



els of refinement (see the discussion in Berger & CoUcla 
1989). To summarize, if shock waves move from high- 



to low-resolution regions, spurious effects associated to 
the mesh structure are undetectable. This is not the case 
in the opposite situation, which causes spurious (though 
small) post-shock oscillations. No te however that for wea k 
shocks the effect is undetectable ( Berger fc Collela 1989|) . 
In cosmology, it is worth mentioning that, since the ba- 
sic features are accretion shocks, we are always in a fa- 
vorable situation: strong shocks originate in high-density 
(high-resolution) regions, and propagates outwards, in a 
low-density (low-resolution) background. To my opinion, 
this fundamental property allows us to use safely adaptive 
mesh technique in cosmological simulations. 



3.6. Spherical Sedov Blast Wave 

We now consider a very difficult test for Cartesian grids 
like the one used in RAMSES: the spherical Sedov test. 
In contrary to the planar, ID case, the spherical blast 
wave is now fully three-dimensional and pretty far from 
the natural geometry of the code. Moreover, as stated be- 
fore, shock waves in 3D are essentially two-dimensional: 



Ax 



dim— 1 



(26) 



where Rs is the curvature radius of the shock. For dim = 3, 
one clearly sees that the number of required cells quickly 
explodes. On the other hand, for the Sedov blast wave 
test, it is more interesting to keep the relative thickness of 
the shock low enough to capture the true solution. As we 
have seen in the last section, if one degrades the resolu- 
tion as the shock propagates outwards, no spurious effects 
are expected. In RAMSES, we can enforce a position- 
dependent spatial resolution by forbidding a given level of 
refinement to be activated if th e ra dius of the cell is larger 
than a given threshold (see § 2.5). The run parameters 
are therefore the followings: the coarse grid size is set to 
= 32 and the maximum level of refinement 
is chosen to be £max = 6. The maximum refinement radius 
for each level is given by 



for < ^ < 5 



(27) 



We use refinement criteria based on pressure gradients 
with Cp = 0.5 and mesh smoothing parameter ricxpand — 
1. The fluid is supposed initially at rest with pq = 1 and 
Pp = 10-5. A total (internal) energy i?o = 1 is deposited in 
the 8 central cells only. The refined mesh is here again built 
before the beginning of the time integration. We assume 
cfl = 0.8 and 7 = 1.4. We use here a single time step for 
all levels, since in this particular case, this is the fastest 
solution. 

Results are shown in Figure |6 and compared to the 
analytical solutions of Sedov (1993 ) for 3 different output 
times. Each quantity represents a volume-average value 
over spherical bins, whose thicknesses correspond to the 
local resolution. Each quantity was also rescaled according 
to the (time-dependent) analytical post-shock values (la- 
beled with an "s") for sake of visibility. Error bars are 
computed using the standard deviation of the numeri- 
cal solution from the mean value in each spherical bin. 
The agreement with the analytical solution is remarkable, 
considering that the mesh has a Cartesian geometry. The 
main departure is found in the velocity profile at a radius 
around 60 % of the shock radius. Similar results were ob- 



tained by Fryxell et al. (2000). An easy way of solving this 
problem would be to lower the pressure gradients thresh- 
old Cp, which would directly increase the resolution in 
this region, at the expense of increasing the total num- 
ber of cells. Oscillations due to spurious mesh refiections 
are not visible in the radial profiles. Direct inspection of 
the 3D data shows that the only systematic effect is the 
departure from spherical symmetry due to the Cartesian 
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Fig. 6. Spherical Sedov Blast Wave Test: Rescaled density, velocity, pressure and volume-averaged refinement level 
as a function of radius for times t = lO""*, 5.7 x 10"'* and 3.2 x 10"'^. Numerical results (solid lines with error bars) 
are compared to the analytical solutions (solid lines). See text for details. 



nature of the mesh. In Figure ^ the volume-averaged re- 
finement level is shown as a function of radius for differ- 
ent times. Due to the maximum refinement radii we used 
(Eq. 11^), the adaptive mesh evolution is also self-similar, 
though in a piecewise constant manner. The total number 
of cells (and therefore the memory used) remains roughly 
constant over the calculation (around 10^ cells, including 
split cells). The interest of using a tree-based approach for 
building the adaptive mesh appears clearly in this test, 
since the mesh structure follows as closely as possible the 
spherical shape of the shock front. 

3.7. Zel'dovich Pancake 

In this section, typical conditions encountered in cosmo- 
logical simulations are addressed using the Zel'dovich pan- 
cake test. This test is widely used to benchmark cosmo- 
logical hydrodynamics codes (Cen 1992; Ryu et al. 1993; 
Bryan & Norman 1997; Teyssier ct al. 199§| ), since it en- 
compasses all the relevant physics (gravity, hydrodynam- 
ics and expansion). It can be thought as a single mode 
analysis of the collapse of random density perturbations, a 
first step towards the study of the fully three-dimensional 
case. The initial conditions are defined for a given start- 
ing redshift Zi in an Einstein-de Sitter universe (fim = 1, 
fib = 0.1), using a sinusoidal density perturbation of unit 
wavelength^ i.e. of the form 



5p 
P 



1 + Ze 
1 



COS (27rx) 



(28) 



where x is the comoving distance to the pancake mid-plane 
(from now on, we always use super-comoving coordinates, 
as defined in § 2_^). The initial velocity field is set accord- 
ing to the linear theory of gravitational instability 



1 



l + ^c 



2^ (1 + Z,)3/2 



sin {2ttx) 



(29) 



The collapse redshift is chosen to be 1 -I- Zc = 10 and the 
initial redshift is 1 + Zi — 100. The initial baryons tem- 
perature was set to a very low arbitrary value, consistent 
with a negligible background temperature. 

The coarse grid is defined by = 32. We use Np = 
256 particles for the dark matter component. The maxi- 
mum level of refinement was set to ^max = 6, correspond- 
ing to a formal resolution of 2048. Two different cases are 
investigated: in the first run, both pressure gradients and 
gas density thresholds (quasi-Lagrangian mesh) are used 
to build the adaptive mesh, with 



Axf, 



and Cf 



0.5 



(30) 



while in the second case, only gas density thresholds are 
used to trigger new refinements. 



Results are shown in Figure |^ for both cases. For the 
first case, the two accretion shock fronts are refined up 
to the maximum refinement level, and are therefore very 
sharp. For the second case, however, the shock fronts are 
not refined at all. The shock waves are traveling outwards, 
from the high-resolution region in the pancake center, to 
the low-resolution background. In light of what have been 
discussed in the previous sections, this explains why no 
oscillations (due to potential spurious reflections at level 
boundaries) are visible. It is worth mentioning that both 
sets of profiles are almost indistinguishable in the center 
of the pancake. This last test is very encouraging, since 
it allows us to avoid refining shock fronts in cosmological 
simulations. The opposite situation would have been dra- 
matic, because of the large filling factor of cosmic shock 
waves (especially in 3D), which would result in a very large 
memory overhead, and because it would trigger coUision- 
ality in the dark matter particles distribution. 

3.8. Spherical Secondary Infall 

The last test, while interesting, is not very stringent, since 
it is very close to the natural, Cartesian geometry of the 
code. An analytical solution describing the fully non-linear 
collapse of spherical density perturbations was found by 



Bertschinger (1985), for both pure dark matter and pure 



baryons fluids. The initial conditions defining the system 
are the foUowings: a completely homogeneous Einstein-de 
Sitter universe (with flm = 1) contains a single mass per- 
turbation SMq at some initial epoch tg. Surrounding this 
initial seed, shells of matter with increasing radius starts 
expanding within the Hubble flow, but finally decouples 
from the expansion at some "turn around" time, and the 
corresponding turn around proper radius, given by 



/ 4 i 

rta[t) = — - 
V Sir ti 



8/9 



Pc 



(31) 



Since in the problem, no other time- or length-scale are 
involved, the overall evolution is self-similar. 



In Kravtsov et al. (1997 ), the secondary infall test was 
successfully passed by the ART code for the pure dark 
matter case. Results obtained by RAMSES are very close 
to the ones obtained by the ART code, which is reassuring, 
since both codes have almost the same N-body solver. 
They are not presented here. Rather, we investigate the 
purely baryonic self-similar infall, so as to validate the 
hydrodynamics solver coupled to gravity and cosmological 
expansion. 

The periodic box is initially filled with a critical density 
cold gas with 7 = 5/3. A single dark matter particle with 
mass rrip = 1/8 is placed as initial seed in the center of the 
computational domain. The coarse grid is defined by = 
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Fig. 7. Zel'dovich Pancake Test: Density, velocity, pressure and refinement level as a function of position from 
the pancake mid-plane at z = 0. The solid line shows AMR results if refinements are activated using both density 
thresholds and pressure gradients. This explains why the 2 accretion shocks are refined. The squares show AMR results 
if refinements are activated using only density thresholds. See text for details. 

Fig. 8. Secondary Infall Test: Rescaled density, velocity, pressure and volume-averaged refinement level as a function 
of radius (in units of coarse cells length) for expansion factors a = Giai, 512ai and 40960.;. Numerical results (solid 
lines with error bars) are compared to the analytical solutions of Bertschinger (1985) (solid lines). See text for details. 



Uy = Uz = 32 and the maximum level of refinement was 
set to fmax = 6, providing us a formal spatial resolution 
of 2048"^. Before the beginning of time integration, the 
mesh is refined around the central seed using a maximum 
refinement radius for each level given by 



Rf. = 2 



3-e 



for < ^ < 5 



(32) 



in units of coarse cell size. The resulting mesh structure 
can be seen in Figure |^ in a volume-averaged radial repre- 
sentation, with roughly 10^ cells in the AMR tree, includ- 
ing split cells. No pressure gradients criterion is used, so 
that shock fronts are not refined. We use a Courant factor 
cfl = 0.5. Starting at expansion factor a; — 10~^, three 
output times were analyzed (a ~ 64ai, 512ai and 40960;). 
The final epoch was reached in 86 coarse time steps only, 
but 5504 time steps at the maximum level of refinement. 

The resulting rescaled density, pressure and velocity 
profiles are plotte d in Figure and c ompared to the ana- 
lytical solution of Bertschinger {198L ) . Error bars are com- 
puted using the standard deviation of the numerical solu- 
tion with respect to the average radial value. The scaling 
relations for velocity, density and pressure are obtained 
using their "turn around" values 



ptait) = (6^Gi2)-l 



Utait) 



rta{t) 



Ptait) = Ptait) 



rtajty^ 



(33) 
(34) 

(35) 



We find a very good agreement between numerical and 
analytical profiles, down to a radius of 2 fine cells, the ac- 
tual resolution limit of the code. As the shock propagates 
outwards, no spurious reflection appears, as expected. 

4. Application to Cosmology: Structure Formation 
in a ACDM universe. 

In this section, results obtained by RAMSES for a N-body 
and hydrodynamical simulation of structure formation in 
a ACDM universe are reported. The box size was set to 
L — 100 h~^ Mpc, as a good compromise between cos- 
mic variance and resolution. The influence of the chosen 
box size on the results are not investigated in this pa- 
per. On the other hand, the convergence properties of the 
solution are examined by varying the mass and spatial 



resolution using 6 different runs, whose parameters are 
listed in Table |l| below. Numerical results are also com- 
pared to analytical results obtained in the framework of 
the halo model. This simple theory predicts various quan- 
tities for both gas and dark matter distributions, and has 
already proven to successfully reproduce results obtained 
in various numerical simulations (see Sect. 4.5). A care- 



ful comparison between the analytical and the numerical 
approach will serve us as a guide to investigate our under- 
standing of structure formation in the universe. 

4.1. Initial conditions 

An initial Gaussian random field was generated for the 
highest resolution run on a 256'^ particle grid and (peri- 
od ic) box len gth L — 100h~^ Mpc. The transfer function 
of Ma (1998| ) for a fiat ACDM universe was used and nor- 
mafized to the COBE data ( [White fc Bunn 1995| ), with 
the following cosmological parameters 



n„, = 0.3 riA = 0.7 ilh = 0.039 



h = 0.7 (Tg = 0.92 



(36) 
(37) 



The high resolution grid was then degraded twice (down 
to 128'^ and 64^) to provide consistent initial conditions 
for our low resolution runs. In this way, a direct com- 
parison between the 3 simulations is made possible. The 
corresponding mass resolution (corresponding to individ- 
ual particle masses for a pure dark matter universe) is 
Mo = 5x10'^ Mq (4x10^° Mq and 3x10" Mq). Particles 
were initially displaced according to the Zel'dovich ap- 
proximation up to a starting redshift 1 + Zi = 72 (51, 
36) for the 256^ (128^, 64^) grid. The initial gas density 
and velocity field was perturbed according to the linear 
theory of gravitational clustering, using the same density 
and displacement fields as for dark matter. The initial gas 
temperature was set to T = 548(1 -f z;)^^ K, in order 
to recover the correct thermal history of baryons after re- 
combination, when neglecting re-ionization. The adiabatic 
index was set to 7 = 5/3, and a fully ionized, primordial 
H and He plasma was considered, with mean molecular 
weight = 0.59. 

4.2. Refinement Strategy 

The main ingredient in the cosmological simulations pre- 
sented here is the refinement strategy. In order to increase 



14 



R. Teyssier: Cosmological Hydrodynamics with Adaptive Mesh Refinement 



Table 1. RAMSES parameters for our six ACDM simu- 
lations. In each case, the box size was set to L = 100 
Mpc. 



Name 


A'^part 


■^max 


A^coll 




Astcp 












h-^kpc 


£ = 


^max 


P064L0 


64^ 





2.6 X 10^ 


1562. 


51 


51 


P128L0 


128^ 





2.1 X 10'^ 


781.3 


107 


107 


P256L0 


256^ 





1.7 X 10'^ 


390.6 


243 


243 


P064L3 


64^ 


3 


5.6 X lO'' 


195.3 


67 


493 


P128L4 


128^ 


4 


5.0 X 10*^ 


48.82 


148 


2259 


P256L5 


256^ 


5 


4.1 X lO'^ 


12.21 


304 


7281 



the spatial resolution within collapsing regions, a quasi- 
Lagrangian mesh evolution was naturally chosen, using 
level dependent dark matt er and gas density thresholds, 
as explained in Section 2_^. To be more specific, the level 
dependent density thresholds for the dark matter compo- 
nent were set to 



Pi = 40- 



nh Ml 







while for the baryonic component, they were set to 
fib Mo 



P£ = 40 



(38) 



(39) 



In this way, a roughly constant number of particles per 
cell (between 5 and 40) is obtained, minimizing both col- 
lisionality and particle discreteness effects. Note that this 
number is higher than the pure dark matter simulations 
performed in Kravtsovct al. (1997D , where 5 particles were 
used to trigger new refinements, instead of 40 in this pa- 
per. As explained above, this choice has basically two rea- 
sons: 1- we prefer to minimize as much as possible the ef- 
fect of particle discreteness effect (Poisson noise) on fluid 
dynamics, 2- since the memory storage is dominated by 
fluid variables, we can increase the number of particles 
per cell by one order of magnitude for a given spatial res- 
olution. The number 40 was finally retained to allow for 
a simple comparison with the more standard refinement 
threshold of Kravtsov et al. (1997 ), namely a factor of 2 



decrease in spatial resolution. 

Shock refinements, as discussed above, was not re- 
tained for these cosmological simulations. This choice has 
two reasons. First, the memory overhead associated to 
shock refinements would have been very large, since shock 
fronts occur everywhere in the hierarchical clustering pic- 
ture. Since shock front are essentially two-dimensional, the 
number of cells required to cover the shock surfaces would 
have been completely out of reach, even for modern com- 
puters. Secondly, refining the mesh in low density regions 
where shock waves eventually propagate would violate the 
non-coUisionality condition for dark matter dynamics. On 
the other hand, the AMR dynamics of shock fronts in this 
case (no shock refinements) was carefully investigated in 
the last sections. It turned out that as soon as shock waves 
travels from high-resolution to low-resolution regions, no 
spurious effects occurs. This last conditions turns out to 



Fig. 9. Top: Energy conservation as a function of expan- 
sion factor for run P256L5 (solid hne), P128L4 (dotted 
line) and P064L3 (dashed line). Bottom: Total number of 
cells in the AMR grid hierarchy (including split cells) di- 
vided by the initial number of coarse cells as a function of 
expansion factor for the same 3 runs. 



be satisfied in cosmological simulations, as discussed in 
the next section. 

The 3 different initial particle grids considered here 
(64^, 128^ and 256^) defines also the coarse level (i = 0) 
of the AMR hierarchy in each run. The maximum level of 
refinement was set to 3, 4 and 5 respectively. This corre- 
sponds to a formal resolution of 512^, 2048^ and 8192^ in 
the highest resolution regions. The corresponding spatial 
resolution is 195, 48 and 12 h^^ kpc. In all cases, adap- 
tive time stepping was activated, with the following time 
step control parameters Ci = C2 = cfl = 0.5. In order 
to measure the advantage of adaptive mesh in cosmologi- 
cal simulations, these 3 runs were also performed without 
refinement (^max — 0). 



4.3. Energy Conservation and Adaptive Mesh 
Evolution 

A standard measure of the quality of a numerical sim- 
ulation is to check for total energy conservation errors. 
Since Euler equations are solved in conservative form in 
RAMSES, the main source of energy conservation er- 
rors comes from the gravitational source terms, for both 
baryons and dark matter components. Figure ^ shows the 
total energy conservation (in the form of the Layzer-Irvine 
conservation equation for an expanding universe) for the 
3 AMR runs. The maximum errors are found to be 2 %, 
1 % and 0.5 % for P064L3, P128L4 and P256L5 respec- 
tively. The maximum error in the energy conservation oc- 
curs when a significant number of refinements are built 
for the first time, around 1 + z 2± 2, 3 and 5 for P064L3, 
P128L4 and P256L5 respectively 

In Figure || is also shown the total number of cells in 
the AMR tree structure (including split cells), in units of 
the number of coarse cells. It is worth noticing that this 
numbers should have remained exactly equal to 1 for a 
strict Lagrangian mesh like the ones described in Gnedin 



1995 ) and Pen (1995 ). At the end of the simulations, as 
clustering develops, the final number of mesh points has 
increased by a factor of 2.5. This overhead is related to 
the mesh smoothing operator (see Sect. |2.1.lD . The mesh 
evolution is therefore only "quasi" Lagrangian. 



4.4. Adaptive Mesh Structure 

The adaptive mesh is dynamically modified at each time 
step during the course of the simulation. Both hydrody- 
namics and N-body solvers take advantage of the increased 
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Fig. 10. Gray scale images of the gas temperature {T > 10^ K), the gas density (p > ilbPc) and the dark matter 
density (p > flmPc) in a planar cut through the computational volume for run P256L5. The mesh structure within 
the plane is plotted in the upper right panel (only octs boundaries are shown for visibility). For sake of comparison, 
the density and temperature maps obtained for run P256L0 (same initial conditions without refinements) are shown 
in the 2 lower panels. 



spatial resolution to improve the accuracy of the solution 
in the refined regions. 

Figure |l^ illustrates this by comparing the gas and 
dark matter density fields in a slice cutting through the 



4.5. The Halo Model 

In order to analyze the results of the simulations in a more 
quantitative way, I will use a powerful analytical theory: 



compu national volume for run P256L5 (256'' particles with 



the so called halo model. Several authors (Seljak 2000; Ma 
fc Fry 2000|; jScoccimarro et al. 200l|; |Cooray et al . 2000l; 



ICooray 200inR'efregier fc Teyssier 2001| ) have recently ex- 
plored the idea that both dark matter and baryons dis- 
tributions can be described by the sum of two contribu- 
tions: (1) a collection of virialized, hydrostatic haloes with 
overdensity > 200 and described by the Press & Schechter 
mass function and (2) a smooth background with overden- 
sity < 10 described by the linear theory of gravitational 
clustering. The purpose of this paper is not to improve 
upon earlier works on this halo model, but rather to use 
it as an analyzing tool for our simulations. Therefore, the 
basic ingredients of the halo model are only briefly recalled 
and will not be discussed in great details. From now on, we 
consider only results obtained at the final redshift z = 0. 
The redshift evolution of the halo model is discussed and 



5 levels of refinements) with the density fields in the same 
slice for run P256L0 (same initial conditions without re- 
finements). Only overdense cells are shown (p > QbPc for 
baryons and p > {Qrn~^b)Pc for dark matter). One clearly 
sees that both gas and dark matter density fields are much 
more dense and clumpy when refinements are activated. 
On the other hand, it is reassuring to see that both sim- 
ulations agree with each other on large scale. 

The corresponding mesh structure is shown in the up- 
per right part of Figure ^ The interest of using a tree- 
based approach for defining the AMR hierarchy is striking: 
the grid structure closely follows the geometry of the den- 
sity field, from a typical filamentary shape at large scale, to 
a more spherical and compact shape in the higher density 
haloes pores. If one examines the central filament conncct - 



ing th e 2 massive haloes in the images of Figure | 10 | , one 
sees that it follows a typical pancake-like structure, with 
2 dark matter caustics on each side of a gas filament. This 
structure, though interesting, is not dense enough to be 
refined by our refinement strategy. We could have lowered 
the density thresholds to trigger new refinements in this 
region, at the price of an increased coUisionality for dark 
matter. This would result in a spurious fragmentation of 
the pancake structure (Melott et al. 1997). 



compared to numerical simulations elsewhere (Refregiei 
et al. 200C|; |Refregier fc Teyssier 200l|). 



Haloes are defined as virialized clump of gas and dark 
matter, with total mass defined as the Virial mass 



vir - —AcPcr: 



3 

vir 



(40) 



The temperature map (T > 10^ K) in the same pla- 
nar cut exhibits the typical flower-like structure of strong 
cosmological accretion shocks around large haloes. These 
strong shocks propagates exclusively in large voids in be- 
tween filaments that intersect each other at the central 
halo position: this is due to the higher gas pressure within 
the filaments that inhibits shock propagation in the direc- 
tion aligned with the filaments. This property of cosmo- 
logical shock waves is of great importance here, since it 
implies that strong shocks propagates almost exclusively 
from high to low resolution regions of the grid. On the 
other hand, weak shocks occurring during sub-halo merg- 
ers along the filaments can enter high resolution regions 
of the mesh as clustering develops. Since the oscillatory 
behavior outline d in Section 3.5 disappears completel y 
for weak shocks ( |Berger fc Collela 1989| ; [Khokhlov 199^) , 
we can conclude that cosmological shocks propagation re- 



where the Virial density contrast Ac is related to the Virial 
overdensity by Ac = rimA. For — 0.3 and z — 0, one 
has A ~ 334 ( |Eke et al. 19981 ). The dark matter follows 
the Navarro, Frenk & White (NFW; [Navarro et al. 1996D 
density profile 



P - 



Ps 



{r/rs){l+r/rsf 



(41) 



whose parameters are connected to the halo Virial mass 

by 



ln(l 



1 



(42) 



mains ree from spurious effccLs assuciaLed Lu Lite adapLive — Schechter 19741) 



The only remaining free parameter (apart from Mvir) is 
the so called concentration parameter c — r^-„/rs. This 
parameter exhibits a good correlation with halo mass in 
numerical simulations that can be fitted analytically to 
match the numerical results ( Ma fc Fry 2000| ). The final 
ingredient is to assume that the halo distribution is de- 
scribed by the Press & Schechter mass function (Press & 



mesh sLrucLure 



The total mass power spectrum can then computed as 
the sum of 2 components P{k) = Pi{k) + P2{k), where 
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Pi(fc) is a non- linear term corresponding to the mass cor- 
relation within halos, and P2{k) is a linear term corre- 
sponding to the mass correlation between 2 halos. Both 
terms have relatively straightforward analytical expres 
sions that are not recalled here (Seljak 2000; Ma & Fry 
2000 |; l^coccimarro et al. 20011 ). 

The gas distribution within halos is supposed to fol- 
low the isothermal hydrostatic equilibrium. The temper- 
ature remains constant within the halo Virial radius, and 
is taken equal to the halo Virial temperature 



(43) 



Note that temperature profiles determined in numerical 
simulations are neither isothermal nor equal to the Virial 
temperature. Therefore, the halo model can only be con- 
sidered as a crude approximation, describing the gas dis- 
tribution in a statistical sense only. Moreover, the tem- 
perature profiles observed in large X-ray clusters is more 
or less affected by cooling flows in the central regions. 
Including all physical ingredients that might affect the 
thermal structure in the core of virialized halos is beyond 
the scope of this paper. Only adiabatic hydrodynamics is 
considered here. 

Solving the hydrostatic equilibrium equation in the 
isothermal case (using the NFW mass distribution) leads 
to the following gas density profile ( Suto et al. 1998| ) 



-b |^ I / \brs/r 



where the dimensionless parameter h is given by 



(44) 



(45) 



The central density po is computed by specifying that the 
total baryons mass within the Virial radius is equal to 

In the next section, the gas pressure power spectrum is 
computed from RAMSES numerical simulations and com- 
pared to the halo model predictions. The pressure power 
spectrum is quite an interesting quantity in cosmology, 
since it is directly related to the Sunyaev-Zeldovich in- 
duced Cosmic Micro wave Background anisotropics angu- 
lar power spectrum ( Sunyaev fc Zeldovich 1980| ; Rephaeli 
1995|}~|lt can be computed within the halo model frame- 
work using the same two terms as for the total mass den- 
sity power spectrum P{k) = Pi{k) + P2{k). Exac t an- 
alytical expressions can be found in pooray (200l| ) and 
Refregier fc Teyssier (200l|). 



4.6. Power Spectra 

In this section, both dark matter density and gas pres- 
sure power spectra are computed and compared to the 
halo model predictions. In order to study the convergence 
properties of the numerical solution, results obtained with 
different mass and spatial resolutions are examined. 



4.6.1. Computing the power spectra 

Computing power spectra for simulations with such high 
dynamical range requires to go beyond traditional meth- 
ods based on regular Cartesian meshes: recall that our 
highest resolution run has a formal resolution of 8192^. We 
use instead a multi-grid method based on a hierarchy of 
nested cubic Ca rtesian grids ( Jenkins ct al. 1998 ; Kravtsov 
|fc Klypin 199£ ) . Each level of the hierarchy corresponds 



to the code AMR levels (from £ = to €max) and cov- 
ers the whole computational volume with cubic grids 
of size n^. At a given level, a dark matter density field 
is computed for each grid using CIC interpolation, and 
all grids are stacked together in a single, co-added density 
field. This density field is then Fourier analyzed using FFT 
technique. From the resulting power spectrum, only modes 
spanning the range 2^ x [fcmin, fcmax] are kept as reliable es- 
timations of the true power spectrum, with fcmin = fcnynq/8 

and /Cinax — ^nynq 

/ 4. The Nyquist frequency fcnynq de- 
pends on the size of the cubic grid, chosen here equal to 
the coarse grid size n^, so that fcnynq = UxTt/L. At the 2 
extreme spatial scales, we have however fcmin = '2,'it/L (for 
^ = 0) and fcmax = Z^^^'fenynq (for I = ^max)- The max- 
imum frequency considered in the present analysis corre- 
sponds therefore to the formal Nyquist frequency of each 
simulation fcmax = Tr/Aa^min (see Table 0). The same pro- 
cedure is applied to the pressure field, except that CIC 
interpolation is no longer needed. 



4.6.2. Results 

The resulting power spectra are shown in Figure ^ for 
the 3 runs with refinements (labeled "AMR") and with- 
out refinements (labeled "PM" ) . For comparison, the dark 
matter and pressure power spectra predicted by the halo 
model are plotted as solid lines. 

The dark matter power spectrum shows a striking 
agreement with the halo model prediction, down to the 
formal resolution limit. Note that the halo model free pa- 
rameters has been tuned in order to recover simulations re- 



sults ( Ma fc Fry 2000 ). Results obtained here are therefore 
consistent with those obtained by other authors (Jenkins 
|ct al. 199S| ; [Kravtsov fc Klypin 1999| ), and can be consid- 
ered as a powerful integrated test of the code. For each 
mass resolution, the power spectrum is plotted up to the 
formal Nyquist frequency. For AMR runs with refinements 
activated, the numerical power spectrum is dominated at 
high wave numbers by the Poisson noise due to particle 
discreteness effects (see the small increase of power around 
fcnynq)- Wc Can couclude that the numerical power spec- 
trum has converged for each run, down to the limit im- 
posed by the finite mass resolution, without being affected 
by the finite spatial resolution. For runs without refine- 
ments, the limited dynamical range has a noticeable effect 
on the resulting power spectrum at much larger scale: the 
spatial resolution is therefore a strong limiting factor in 
this case. 
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Fig. 11. Dark matter density (left panels) and gas pressure (right panels) power spectra for RAMSES runs with 
refinements (upper panels) and without refinements (lower panels). In each plot, the solid line corresponds to the halo 
model prediction, while the dotted line (dashed and dot-dashed) corresponds to numerical results with 256'^ (128"^ and 
64^) particles. Label "AMR" stands for AMR runs (P256L5, P128L4 and P64L3), while label "PM" stands for runs 
without refinement (P256L0, P128L0 and P64L0). For each power spectrum, the curve ends at the Nyquist frequency 
corresponding to the formal resolution of the simulation (Axmin in Table 



Let us now examine the gas pressure power spectrum 
(Fig. ^l]). The overall agreement of the numerical results 
with the halo model predictions is relatively good: the cor- 
rect behavior is captured at all scales within a factor of 
2 for our highest resolution run (P256L5). Note that the 
halo model has no free parameters for the gas distribu- 
tion, as soon as the dark matter parameters are held fixed. 
At large scale, the numerical power spectrum appear to 
converge to the halo model predictions (run P128L4 and 
run P256L5 give exactly the same results). Note that a 
rather high mass resolution is needed for this convergence 
to occur (> 128'^ particles), as opposed to the dark matter 
density power spectrum for which the correct large scale 
power is recovered even with 64'^ particles. At intermedi- 
ate scales (around 1 h Mpc~^), the halo model predictions 
slightly underestimate the pressure power spectrum. Since 
numerical results have also converged at these scales, this 
discrepancy might be due to the fact that intermediate 



Table 2. Global properties of the 5 largest halos ex- 
tracted from the highest resolution run P256L5. 



Name 






c 


/3flt 


^corc 




h-^M© 


h-"Mpc 






h-^kpc 


Cluster 1 


6.97 X 10"* 


1.82 


9.5 


0.85 


144.2 


Cluster 2 


5.09 X 10" 


1.64 


5.9 


0.80 


215.5 


Cluster 3 


5.07 X 10" 


1.63 


7.3 


0.80 


202.3 


Cluster 4 


4.52 X 10" 


1.57 


4.9 


0.64 


147.0 


Cluster 5 


4.29 X lO" 


1.55 


9.4 


0.84 


159.7 



possible thanks to the large dynamical range obtained in 
our simulation (Axmin = 12 h^^ kpc), although it was not 
optimized to study individual halo properties. The next 
step to go beyond what is presented here would be to 
perform so called "zoom simulations" , with nested, higher 
resolution, initial conditions particle grids centered on sin- 



densitj regions (10 < p < 200) are completely neglected 
in the halo model. These regions are believed to be com- 
posed of warm filaments, whose pressure obviously cannot 
be completely neglected. 

At small scales, the situation remains quite unclear. 
On one hand, one clearly sees on Figure |l^ that an in- 
creased dynamical range has a dramatic effect on the re- 
sulting pressure power spectra. The power is much higher 
on small scales for AMR runs than for runs without re- 
finements. On the other hand, for AMR runs, the con- 
vergence of the numerical results to the "true" solution is 
not as fast as that of the dark matter power spectrum. 
Some hints of convergence between run P128L4 and run 
P256L5 can be seen on Figure [l^. Indeed, without refine- 
ments (runs P64L0, P128L0 and P256L0), the cut off in 
the pressure power spectrum is directly proportional to 
the spatial resolution of the simulation. The same is true 
between runs P64L3 and P128L4, while for run P256L5, 
the effect of spatial resolution appears to weaken slightly. 
More interesting is the large discrepancy in slope at large k 
between the halo model and the solution obtained by run 
P256L5. As discussed in the next section, this is probably 
due to the assumption of isothermality in the halo model, 
which is ruled out by simulation results for individual halo 
temperature profiles. 

4. 7. Individual Halos Structure 

The internal structure of the 5 largest halos found in the 
highest resolution run (P256L5) is now examined in great 
detail. It is worth mentioning that this analysis is made 



gle halos (Bryan fc Norman 1997 ; Eke et al. 1998; Yoshida 
et al. 20001 ; [Abel et al. 2000| ). The main advantage of the 



present "brute force" approach is to combine both large 
and small scale results in the analysis. In order to study 
the convergence properties of individual halo profiles, re- 
sults obtained for runs P256L5, P128L4 and P64L3 are 
compared. Recall that a direct comparison of the same 
halo at different mass and spatial resolution is possible, 
since the same initial conditions were used (and degraded 
to the correct mass resolution) for each run. Runs without 
refinement (P256L0, P128L0 and P64L0) are discarded 
from this analysis, because they completely lack the nec- 
essary dynamical range to resolve the internal structure 
of individual halos. 

Haloes were detected in the dark matter particles dis- 
tribution at the final output time {z — 0) using the 



Spherical Overdensity algorithm (Lacey & Cole 1993), 
with overdensity threshold A = 334. Only the 5 most 
massive haloes of the resulting mass function are consid- 
ered for the present analysis. Their global properties are 
listed on Table |^. For each halo, the center is defined as 
the location of the maximum in the dark matter density 
field. For regular, relaxed haloes, this definition also corre- 
sponds to the maximum in the baryons density field, and 
to the halo center of mass. This is however not the case for 
irregular, not yet relaxed haloes, for which this definition 
of the halo center is less robust. 

Cubic regions 6.25 h~^ Mpc aside are then ex- 
tracted around each halo center. The resulting projected 
color maps for various relevant quantities are shown in 
Figure |l^. Clusters 1 and 5 appear to be the most relaxed 
halos of our sample, while clusters 2, 3 and 4 show more 
substructures and irregularities within their Virial radii. 
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Fig. 12. Color maps showing various projected quantities for the 5 most massive halos extracted from run P256L5. 
The projected volume is in each case a cube of 6.25 h^^Mpc aside. The color coding is based on a logarithmic scale 
for each plot. From top to bottom: 1- Projected dark matter particle distribution (the color coding corresponds to 
the local particle density). 2- X-ray emissivity map. 3- X-ray (emission weighted) temperature map. 4- Sunyaev- 
Zeldovich decrement parameter (equivalent to the integrated pressure along the line of sight). 5- Projected adaptive 
mesh structure (only oct boundaries are shown). 



Fig. 13. Radial profiles for the 5 most massive halos showing various quantities averaged along radial bins. In each 
plot, the dotted line (dashed and dot-dashed) comes from run P256L5 (P128L4 and P64L3). From top to bottom: 1- 
Dark matter overdensity profile (the best-fit NFW analytical profile is shown as a solid line). 2- Gas overdensity profile 
(the corresponding "hydrostatic isothermal model" analytical profile is also shown as a solid line). 3- Mass averaged 
temperature profile (the corresponding "hydrostatic beta model" analytical profile is also shown as a solid line). 4- 
Pressure profile (the corresponding "hydrostatic isothermal model" analytical profile is also shown as a solid line). 5- 
Volume averaged refinement levels. 



The adaptive mesh structure, also shown in Figure |l^, 
closely matches the clumpy structure of each halo. Note 
that the maximum level of refinement (£ — 5) is acti- 
vated in the halo cores, where the formal spatial resolu- 
tion reaches 12 h~^ kpc (barely visible in Fig. ^2|). The 
physical properties of these 5 haloes are now discussed 
quantitatively. 

4.7.1. Dark Matter Distribution 

The projected dark matter particles distribution is shown 
in Figure |l^, with a color coding corresponding to the lo- 
cal particle overdensity. The dark matter distribution is 
far from being smooth and spherically symmetric. It is 
however interesting to compute the radial density profile 
and compare it to the NFW analytical prediction. The 
result is shown in Figure |l^, for our 3 different mass res- 
olutions (runs P256L5, P128L4 and P64L3). The density 
profile obtained in the highest resolution run (P256L5) 
was fitted to the NFW profile, using the concentration 
parameter c as fitting parameter. Best fit values are listed 
in Table ^: they are full y consistent with expected values 
for a ACDM universe ( Kravtsov et al. 1997 ; Eke et al 
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The quality of the fit is impressive for clusters 1 and 
5, which are also the more relaxed haloes of our sample. 
The poorest fit was obtained for cluster 4, with deviations 
as large as 50 % at a radius of 150 h~^ kpc. A close ex- 
amination of the corresponding map in Figure |l^ confirms 
that this halo is poorly relaxed in its central region. By 
comparing the profiles obtained for different mass resolu- 
tion, one sees that the numerical profiles agree with each 



other down to their resolution limit. This is in complete 
agreement with our conclusion concerning the dark matter 
power spectrum: simulated power spectra match closely 
the halo model prediction down to their formal Nyquist 
frequency. 

4.7.2. Baryons Distributions 

The baryons distribution within the selected haloes is sim- 
ilar on large scales to the dark matter distribution. In 
Figure |l^, simulated X-ray emission maps (using Lx oc 
rig) are good tracers of the gas overdensity projected along 
the line of sight. One can notice however that the hot gas is 
more smoothly distributed than dark matter. This is even 
more striking in the central region of all clusters, where 
the gas density reaches a plateau, reminiscent of a /3 model 
density profile. It is worth noticing that overdensities (sub- 
structures) in the gas distribution usually appears as cold 
spots in the X-ray (emission weighted) temperature map. 
On the other hand, the cluster cores are significantly hot- 
ter than the surrounding gas in most cases. We will come 
back to this point later. 

Using the halo center as defined above, gas overden- 
sity, mass weighted temperature and pressure profiles were 
computed as a function of radius and plotted in Figure ^ 
For that purpose, conserved variables such as mass and 
internal energy are averaged into radial bins of increas- 
ing thickness, starting from the formal resolution (12 h~^ 
kpc for run P256L5) up 3 times this value at the Virial 
radius. The volume averaged refinement level is also plot- 
ted in Figure |l^, giving some hints of the effective spa- 
tial resolution as a function of radius. Based on the re- 
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suits obtained during the Spherical Secondary Infall test 
presented in 



3.8 



the actual resolution of the code cor- 
responds roughly to twice its formal resolution. For run 
P256L5 (P128L4 and P64L3), this gives a Umiting radius 
of 24 (96 and 384) h^^ kpc, above which numerical results 
are fully reliable. 

Since dark matter density profile are well fitted for each 
halo by the NFW analytical profile, the corresponding 
gas density profile can be computed using Equation (Q). 
Recall that in the halo model, the gas temperature is as- 
sumed to remain constant within the Virial radius, and 
equal to the Virial temperature (Eq. ^). This is obvi- 
ously not the case for our simulated clusters (see Fig. |l3|) , 
and explains why the halo model profile is much more 
peaked in the central region than for simulated profiles. 
For the pressure profiles, the situation is however less dra- 
matic, though still unsatisfactory. It is interesting to no- 
tice that here again, the behavior of the pressure profiles is 
fully consistent with previous results concerning the pres- 



Since both gas and dark matter density profiles are 
now determined to a good accuracy, it is possible to per- 
form a consistency check and compute the temperature 
profile resulting from the assumption of hydrostatic equi- 
librium. For that purpose, the analytical framework de- 
veloped by Varnieres & Arnaud (2001) is exactly what 
is needed here: assuming a NFW profile for dark matter 
and a (3 model for baryons, they derived analytically the 
corresponding hydrostatic temperature profile. This tem- 
perature profile is shown for each cluster in Figure ^ 
The agreement with numerical results is good: temper- 
ature rising towards the halo center is therefore a direct 
consequence of hydrostatic equilibrium. After closer exam- 
ination, the small (20 %) disagreement observed in clusters 
2, 3 and 4 is due to the poorer fit of the NFW formula to 
the dark matter simulated density profiles. 

This non constant behavior of the halo temperature 
profile have been already noticed by several authors in 
nume rical simulations ( Bryan & Norman 1997]; |Eke et al 



sure power spectrum: the isothermal halo model overes- 
timates the pressure power on small scales. This trans- 
lates in a much steeper slope at low radii, as for the pres- 
sure power spectrum at large k. Moreover, numerical re- 
sults for the gas distribution within haloes have clearly 
not converged yet, although the dependance of the com- 
puted profiles with respect to the spatial resolution seems 
to weaken slightly between run P128L4 and run P256L5, 
in exactly the same way as the pressure power spectrum 
did in S 4.6.2. Conclusions are therefore similar: numerical 



1998| ; Loken et al. 2000 ) and also in X-ray obser vations 
of large galaxy clusters ( Markevitch et al. 1998| ). Note 



results show good evidence of converging at scales greater 



than 53 h ^ kpc. Similar conclusions were obtained by 



however that more physics need to be included in current 
numerical simulations before performing a reliable com- 
parison to observations. On the other hand, the idealized 
case of adiabatic gas dynamics is still of great theoretical 
interest: one can hope to find a self-consistent description 
of the gas distribution within haloes using such an ap- 
proach. The main ingredient to extend the halo model in 
this framework is to determine the typical core radius (as 
well as /3fit) as a function of halo mass and redshift. Eka 



Bryan & Norman (1997) for a "zoom" simulation of a sin- 
gle cluster. These authors obtained very similar gas den- 
sity and temperature profiles, with quite the same conver- 
gence properties as the one obtained here. Consequently, 
the isothermal halo model is likely to fail at scales less 
than 250 h'^ kpc. 

4.7.3. Beyond the Isothermal Halo Model ? 



As noted by several authors (Bryan & Norman 1997; Eke 



et al. 1998), the typical gas density profile obtained in 
numerical simulations is much more accurately described 
by the /3 mod el analytical solution ( [Cavaliere fc Fus^ 
Femiano 1976 ) 

p = po[l + (r-/reo..e)^]"'*"^' (46) 

Although the numerical results presented here have not 
fully converged yet, it is worth exploring an alternative 
to the isothermal halo model. The gas overdensity profile 
obtained in run P256L5 was thus fitted with the /? model 
formula, using both rcoic and /3fit as fitting parameters. 
Best fit values are listed in Table ^ and are consis tent 



with typica l numbers quoted by other authors (e.g. |Eke 
et al. 199^ ). It is worth mentioning that the quality of 



the fit is excellent in each case, except for cluster 4, as 
expected from the previous analysis on the dark matter 
distribution. 



et al. (1998) have already initiated this challenging task 
using high resolution (zoom) SPH simulations for large 
mass haloes (M ~ 10^^ h^^ Mq) , but more work need to be 
done to probe a larg er mass and redshift range. Using the 
analytical formula of Varnieres fc Arnaud (2001 ), it would 
be straightforward to determine the corresponding tem- 
perature profile, and thus to complete the description of 
the baryons component within this extended halo model. 
This ambitious project will not be addressed here, but will 
be considered in a near future. 

5. Conclusions and Future Projects 

A new N-body/hydrodynamical code, RAMSES, has been 
presented and tested in various configurations. RAMSES 
has been written in FORTRAN 90 and optimized on a 
vectorized hardware, namely the Fujitsu VPP 5000 at 
CEA Grenoble. A parallel version was also impemented 
on shared-memory systems using OpenMP directives. A 
distibuted memory version of RAMSES is currently un- 
der construction using a domain decomposition approach. 

The main features of the RAMSES code are the fol- 
lowings: 

1- the AMR grid is built on a tree structure, with new 
refinements dynamically created (or destroyed) on a cell- 
by-cell basis. This allows greater flexibility to match com- 
plicated flow geometries. This property appears to be es- 
pecially relevant to cosmological simulations, since clumpy 
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structures form and collapse everywhere within the hier- 
archical clustering scenario. 

2- the hydrodynamical solver is based on a second or- 
der Godunov method, a modern shock capturing method 



that ensures exact total energy conservation, as soon as 
gravity is not included. Moreover, shock capturing relies 
on a Riemann solver, without any artificial viscosity. 

3- the refinement strategy that was retained for cosmo- 
logical simulations is based on a "quasi-Lagrangian" mesh 
evolution. In this way, the number of dark matter particles 
per cell remains roughly constant, minimizing two-body 
relaxation and Poisson noise. On the other hand, this re- 
finement strategy is not optimal for baryons, since one 
neglects to refine shock fronts (this would have been too 
costly anyway). It has been carefully shown that in this 
case, as soon as strong shocks propagate from high to low 
resolution regions of the grid, no spurious effects appear. 

The code has been tested in standard gas dynamical 
test cases (Sod's test and Sedov's test), but also for in- 
tegrated cosmological tests, like Zel'dovich pancake col- 
lapse or Bertschinger spherical secondary infall. It has 
been shown that the actual resolution limit of the code 
is equal to roughly twice the cell size of the maximum 
refinement level. 

The RAMSES code has been finally used to study the 
formation of structures in a low-density ACDM universe. 
A careful convergence analysis has been performed, using 
the same initial conditions with various mass and spatial 
resolutions, for a fixed box size L = 100 h~^ Mpc. The 
initial number of cell (at the coarse level) was set equal to 
the initial particle grid (64^, 128^ or 256^), for a final num- 
ber of cells only 2.5 larger. The formal spatial resolution 
in the largest run was 8192^ or 12 h^^ kpc comoving. 

Numerical results have been compared to the analyti- 
cal predictions of the so called halo model, for both dark 
matter and gas pressure power spectra, as well as indi- 
vidual haloes internal structure. A good agreement was 
found between the halo model and the numerical results 
for dark matter, down to the formal resolution limit. For 
the baryons distribution, numerical results show some ev- 
idence of converging at scales greater than 50 h~^ kpc 
for our highest resolution run. The halo model reproduces 
simulations results only approximatively (within a factor 
of 2) at these scales. 

A simple extension of the halo model for the fluid com- 
ponent has been proposed. The idea is to assume that the 
average gas density profile within haloes is described by 
a /3 model, whose parameters still need to be determined 
from first principles or from numerical simulations and for 
a rather large mass range, which is far beyond the scope 
of this paper (see however Eke et al. 1998). It is then pos- 
sible to deduce from hydrostatic equilibrium an analytical 



as initial conditions, in order to improve mass and spatial 
resolutions inside individual haloes. This approach seems 
indeed very natural within the AMR framework, and has 
already proven to be successful in recent attempts (Bryan 
fc Noi-man 1997|; [Abel et al. 2000|). Future efforts in the 



temperature profile (Varnieres & Arnaud 2001) that ac 



curately matches the simulation results presented in this 
paper, and should therefore improve considerably the halo 
model. 

Extending the current work to "zoom" simulations is 
currently under investigation, using a set of nested grids 



RAMSES code development will be however more focused 
on including more physics in the description of the gaseous 
component, like cooling, star formation and supernovae 
feedback. 
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